/***************************************************************************
 *   Copyright (C) 2006-2011 by  Institute of Combustion Technology        *
 *   d.mayer@itv.rwth-aachen.de                                            *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/
#include <iostream>

#include <stdio.h>
#include <math.h>

#include "solve.h"
#include "data.h"

//------------------------------------------------------
bool solve(sData* data)
{
   std::cout << "\nCalculation:\n------------\n";

   // FIXME
   
   calcFlux(data);
   
   //was ist Abbruchbedingung
   //was ist die zu lösende GLG
   

   std::cout << "\n";
   return true;
}

//------------------------------------------------------
void calcFlux(sData* data)
{
   static sFace* curFace=0;

   // compute numerical flux of each face
   for(int fId=0; fId<data->nFaces; fId++) {
      curFace=&data->faces[fId];
      // FIXME
	  double u=curFace->u;
	  double v=curFace->v;
	  double rho = data->rho;
	  double lambda = data->alpha;
	  
	  if(curFace->dx == 0)
	  {
		  double phi;
		  if(u<0){phi=curFace->neighCells[M]->phi;}else{phi=curFace->neighCells[P]->phi;}
		  double dPhiDy = (curFace->neighCells[M] - curFace->neighCells[P]) / (2 * curFace->dy);
		  
		  curFace->numFlux[0]=rho*u*phi-lambda*dPhiDy;
		  curFace->numFlux[1]=0;
	  }
	  
	  if(curFace->dy == 0)
	  {
		  double phi;
		  if(v<0){phi=curFace->neighCells[M]->phi;}else{phi=curFace->neighCells[P]->phi;}
		  double dPhiDx = (curFace->neighCells[M] - curFace->neighCells[P]) / (2 * curFace->dx);
		  
		  curFace->numFlux[0]=0;
		  curFace->numFlux[1]=rho*u*phi-lambda*dPhiDx;
	  }
	  
   }
}

